function r = gauss_int2d_tri(function_f, V)
gauss_points_tri = [0.11270166537925829786104259255808, 0.1;
    0.11270166537925829786104259255808, 0.44364916731037085106947870372096;
    0.11270166537925829786104259255808, 0.78729833462074172434341789994505;
    0.5, 0.056350832689629148930521296279039;
    0.5, 0.25;
    0.5, 0.44364916731037085106947870372096;
    0.88729833462074170213895740744192, 0.01270166537925830792243875322356;
    0.88729833462074170213895740744192, 0.056350832689629148930521296279039;
    0.88729833462074170213895740744192, 0.1];
gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
J_det = (V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1));
N = reference_basis_2d("P1");
r = 0;
for i = 1:9
    xi = gauss_points_tri(i,1);
    eta = gauss_points_tri(i,2);
    X = N{1}(xi,eta)*V(1,:) + N{2}(xi,eta)*V(2,:) + N{3}(xi,eta)*V(3,:);
    r = r + gauss_weights_tri(i)*function_f(X(1), X(2));
end
r = r * abs(J_det);
end